#include "rsagricuturedroughtwaterlogging.h"
#include <QString>
#include <QObject>
#include <QThread>

RSAgricutureDroughtWaterlogging::RSAgricutureDroughtWaterlogging()
{
    m_ImgWidth = 0;
    m_ImgHeight = 0;

    m_pNDVI = nullptr;
    m_pVSWI = nullptr;
    m_pFloodArea = nullptr;

    m_NewInputImg = "";
    m_Flag = -1;
}

RSAgricutureDroughtWaterlogging::~RSAgricutureDroughtWaterlogging()
{
    if ( nullptr != m_pNDVI )
    {
        delete []m_pNDVI;
        m_pNDVI = nullptr;
    }
    if ( nullptr != m_pVSWI )
    {
        delete []m_pVSWI;
        m_pVSWI = nullptr;
    }
    if ( nullptr != m_pFloodArea )
    {
        delete []m_pFloodArea;
        m_pFloodArea = nullptr;
    }
}

bool RSAgricutureDroughtWaterlogging::setNewInputImg(const char *pNewInputImg)
{
    if ( nullptr == pNewInputImg )
    {
        QString strMsg = QObject::tr("string is NULL");
        m_Msg = strMsg.toUtf8().data();
        return false;
    }

    m_NewInputImg = pNewInputImg;

    return true;
}

void RSAgricutureDroughtWaterlogging::setFlag(const int nFlag)
{
    m_Flag = nFlag;
}

void RSAgricutureDroughtWaterlogging::setOutputImgFormat(const char *pImgFormat)
{
    m_pszFileout = pImgFormat;
}

double RSAgricutureDroughtWaterlogging::getNDVI() const
{
    if ( nullptr == m_pNDVI )
    {
        return 0.0;
    }

    return *m_pNDVI;
}

bool RSAgricutureDroughtWaterlogging::disPose()
{
    //开始
    QString strMsg = QObject::tr("start");
    setDisPoseState( 0 , strMsg.toUtf8().data());

    for ( int i = 0 ; i < 100 ; i++ )
    {
        QThread::msleep( 1 );
        strMsg = QObject::tr("NO:%1").arg(i + 1);
        setDisPoseState( i , strMsg.toUtf8().data());
    }

    //注册文件格式
    GDALAllRegister();
    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");

    //数据集
    GDALDataset *poDatasetold;
    GDALDataset *poDatasetnew;

    //使用只读方式打开第一期图像
    poDatasetold = (GDALDataset*) GDALOpen( m_InputImgVec[0].c_str(),GA_ReadOnly );
    if( NULL == poDatasetold )
    {
        //无法打开第一个文件，返回-1
        strMsg = QObject::tr("open input image fail");
        setDisPoseState( 0 , strMsg.toUtf8().data());
        m_Msg = strMsg.toUtf8().data();
        return false;
    }

    int rasterband = poDatasetold->GetRasterCount();
    if (0 == m_Flag)
    {
        if (2 != rasterband)
        {
            //关闭文件
            GDALClose((GDALDatasetH)poDatasetold);

            //第一个文件波段较少，返回-2
            strMsg = QObject::tr("the bands number of input image ");
            setDisPoseState( 0 , strMsg.toUtf8().data());
            return false;

        }
    }
    else
        if (1 == m_Flag)
        {
            if (5 != rasterband)
            {
                //关闭文件
                GDALClose((GDALDatasetH)poDatasetold);

                //第一个文件波段较少，返回-2
                strMsg = QObject::tr("the bands number of input image ");
                setDisPoseState( 0 , strMsg.toUtf8().data());
                return false;

            }
        }
        else
            if(4 != rasterband)
            {
                //关闭文件
                GDALClose((GDALDatasetH)poDatasetold);

                //第一个文件波段较少，返回-2
                strMsg = QObject::tr("the bands number of input image ");
                setDisPoseState( 0 , strMsg.toUtf8().data());
                return false;
            }

    //使用只读方式打开第二期图像
    poDatasetnew = (GDALDataset*) GDALOpen( m_NewInputImg.c_str(),GA_ReadOnly );
    if( NULL == poDatasetnew )
    {
        //关闭文件
        GDALClose((GDALDatasetH)poDatasetold);

        //打开第二个文件失败，返回-3
        strMsg = QObject::tr("the second input image open fail ");
        setDisPoseState( 0 , strMsg.toUtf8().data());
        return false;
    }

    rasterband = poDatasetnew->GetRasterCount();
    if (0 == m_Flag)
    {
        if (1 != rasterband)
        {
            //关闭文件
            GDALClose((GDALDatasetH)poDatasetold);
            GDALClose((GDALDatasetH)poDatasetnew);

            //第一个文件波段较少，返回-2
            strMsg = QObject::tr("the bands number of input image ");
            setDisPoseState( 0 , strMsg.toUtf8().data());
            return false;

        }
    }
    else
        if (1 == m_Flag)
        {
            if (5 != rasterband)
            {
                //关闭文件
                GDALClose((GDALDatasetH)poDatasetold);
                GDALClose((GDALDatasetH)poDatasetnew);

                //第一个文件波段较少，返回-2
                strMsg = QObject::tr("the bands number of input image ");
                setDisPoseState( 0 , strMsg.toUtf8().data());
                return false;

            }
        }
        else
            if(4 != rasterband)
            {
                //关闭文件
                GDALClose((GDALDatasetH)poDatasetold);
                GDALClose((GDALDatasetH)poDatasetnew);

                //第二个文件波段数较少，返回-4
                strMsg = QObject::tr("the bands number of input image ");
                setDisPoseState( 0 , strMsg.toUtf8().data());
                return false;
            }

    //图像的大小
    int nXsizeold = poDatasetold->GetRasterXSize();
    int nYsizeold = poDatasetold->GetRasterYSize();
    int nXsizenew = poDatasetnew->GetRasterXSize();
    int nYsizenew = poDatasetnew->GetRasterYSize();

    //获取第一个图像的投影及坐标转换参数
    double GTold[6];
    poDatasetold->GetGeoTransform(GTold);
    const char * projold = poDatasetold->GetProjectionRef();

    //获取第二个图像的投影及坐标转换参数
    double GTnew[6];
    poDatasetnew->GetGeoTransform(GTnew);
    const char * projnew = poDatasetnew->GetProjectionRef();

    //判断两幅图像的投影是否一致
    if (0 < m_Flag)
    {
        if (0 != strcmp(projold,projnew))
        {
            //关闭文件
            GDALClose((GDALDatasetH)poDatasetold);
            GDALClose((GDALDatasetH)poDatasetnew);

            //投影不一致，返回-5
            strMsg = QObject::tr("the two of image project  ");
            setDisPoseState( 0 , strMsg.toUtf8().data());
            return false;
        }
    }
    std::string pro = projold;
    int npro = pro.size() + 1;
    //投影信息
    char * proj = new char[npro];
    strcpy(proj,projold);
    //地理变换矩阵
    double *GT = new double[6];

    //判断两幅图像是否有公共区域
    if (GTold[0] >= GTnew[0])
        GT[0] = GTold[0];
    else
        GT[0] = GTnew[0];
    if(GTold[3] <= GTnew[3])
        GT[3] = GTold[3];
    else
        GT[3] = GTnew[3];

    GT[1] = GTold[1];
    GT[2] = GTold[2];
    GT[4] = GTold[4];
    GT[5] = GTold[5];

    double xold = GTold[0] + nXsizeold * GTold[1];
    double yold = GTold[3] + nYsizeold * GTold[5];

    double xnew = GTnew[0] + nXsizenew * GTnew[1];
    double ynew = GTnew[3] + nYsizenew * GTnew[5];

    double x,y;
    if (xold <= xnew)
        x = xold;
    else
        x = xnew;
    if (yold >= ynew)
        y = yold;
    else
        y = ynew;

    m_ImgWidth = (x - GT[0])/GT[1];
    m_ImgHeight = (y - GT[3])/GT[5];

    if (m_ImgWidth <= 0 || m_ImgHeight <= 0)
    {
        //关闭文件
        GDALClose((GDALDatasetH)poDatasetold);
        GDALClose((GDALDatasetH)poDatasetnew);

        //没有共同区域，返回-6
        strMsg = QObject::tr("the bands number of input image ");
        setDisPoseState( 0 , strMsg.toUtf8().data());
        return false;
    }
    //计算旱情
    if (0 == m_Flag)
    {
        //计算NDVI
        m_pNDVI = calculateNDVI(poDatasetold,nXsizeold,nYsizeold);

        //申请存放VSWI内存
        m_pVSWI = (unsigned char*) CPLMalloc(sizeof(unsigned char)*m_ImgWidth*m_ImgHeight);

        //中间变量
        float *temp = (float*) CPLMalloc(sizeof(float)*m_ImgWidth*m_ImgHeight);

        //读取地表温度数据
        GDALRasterBand *poBand = poDatasetnew->GetRasterBand(1);
        unsigned short *LST = (unsigned short*) CPLMalloc(sizeof(unsigned short)*nXsizenew*nYsizenew);

        if (NULL == m_pNDVI || NULL == m_pVSWI || NULL == LST)
        {
            return false;
        }

        poBand->RasterIO(GF_Read, 0, 0, nXsizenew,nYsizenew,LST, nXsizenew,nYsizenew, GDT_UInt16, 0, 0 );

        //温度系数
        float nscale = 0.02;

        //遍历
        strMsg = QObject::tr("Calculate ");
        for (int i = 0; i < m_ImgHeight; i++)
        {
            //按行读取图像数据

            for (int j = 0; j <m_ImgWidth; j++)
            {
                //当前点在NDVI中的位置
                int jold = (GT[0] - GTold[0] + j * GT[1])/GTold[1];
                int iold = (GT[3] - GTold[3] + i * GT[5])/GTold[5];

                //当前点在LST中的位置
                int jnew = (GT[0] - GTnew[0] + j * GT[1])/GTnew[1];
                int inew = (GT[3] - GTnew[3] + i * GT[5])/GTnew[5];

                //计算vswi
                if (LST[inew * nXsizenew + jnew] == 0)
                    temp[i * m_ImgWidth + j] = 0;
                else
                    temp[i * m_ImgWidth + j] = m_pNDVI[iold * nXsizeold + jold]/(nscale * LST[inew * nXsizenew + jnew]);
            }

            setDisPoseState( (i + 1) * 1000 / m_ImgHeight , strMsg.toUtf8().data());
        }

        //测试进程条
        for ( int i = 0 ; i < 1000 ; i++ )
        {
            setDisPoseState( i , strMsg.toUtf8().data());
            QThread::sleep(500);
        }

        //计算最大最小值
        float nmax = getMax(temp);
        float nmin = getMin(temp);

        for (int i = 0;i < m_ImgWidth * m_ImgHeight; i++)
            m_pVSWI[i] = (unsigned char)((temp[i]-nmin)/(nmax-nmin)*256);

        //存储结果图像
        saveImageproj(GT,proj,m_pVSWI,m_pszFileout.c_str());

        //释放内存
        CPLFree(temp);
        CPLFree(LST);
    }
    //计算水淹区...
    if (m_Flag > 0)
    {
        bool *floodold,*floodnew;
        m_pFloodArea = (unsigned char*) CPLMalloc(sizeof(unsigned char)*m_ImgWidth*m_ImgHeight);
        if (1 == m_Flag)
        {
            //计算一期水域
            floodold = modisFlood(poDatasetold,nXsizeold,nYsizeold);

            //计算二期水域
            floodnew = modisFlood(poDatasetnew,nXsizenew,nYsizenew);
        }
        else
            if (2 == m_Flag)
            {
                //计算一期水域
                floodold = landsat8Flood(poDatasetold,nXsizeold,nYsizeold);

                //计算二期水域
                floodnew = landsat8Flood(poDatasetnew,nXsizenew,nYsizenew);

            }

        strMsg = QObject::tr("Calculate ");
        for (int i = 0; i < m_ImgHeight; i++)
        {
            for (int j = 0; j <m_ImgWidth; j++)
            {
                int jold = (GT[0] - GTold[0] + j * GT[1])/GTold[1];
                int iold = (GT[3] - GTold[3] + i * GT[5])/GTold[5];

                int jnew = (GT[0] - GTnew[0] + j * GT[1])/GTnew[1];
                int inew = (GT[3] - GTnew[3] + i * GT[5])/GTnew[5];

                if (floodold[iold * nXsizeold + jold] == 0 && floodnew[inew * nXsizenew + jnew] == 1)
                    m_pFloodArea[i * m_ImgWidth + j] = 255;
                else
                    m_pFloodArea[i * m_ImgWidth + j] = 0;
            }
            setDisPoseState( (i + 1) * 1000 / m_ImgHeight , strMsg.toUtf8().data());
        }

        //关闭文件
        if (NULL != poDatasetold)
            GDALClose((GDALDatasetH)poDatasetold);
        if (NULL != poDatasetnew)
            GDALClose((GDALDatasetH)poDatasetnew);

        //释放内存
        if (NULL != floodold)
        {
            delete []floodold;
            floodold = NULL;
        }
        if (NULL != floodnew)
        {
            delete []floodnew;
            floodnew = NULL;
        }
        //结果写入文件
        if (-1 == saveImageproj(GT,proj,m_pFloodArea,m_pszFileout.c_str()))
        {
            if (NULL != GT)
            {
                delete []GT;
                GT = NULL;
            }
            if (NULL != proj)
            {
                delete []proj;
                proj = NULL;
            }
            //写文件失败，格式不正确
            return false;
        }
    }



    ////释放内存
    if (NULL != GT)
    {
        delete []GT;
        GT = NULL;
    }
    if (NULL != proj)
    {
        delete []proj;
        proj = NULL;
    }

    strMsg = QObject::tr("end");
    setDisPoseState( 1000 , strMsg.toUtf8().data());

    return true;
}

float *RSAgricutureDroughtWaterlogging::calculateNDVI(GDALDataset *poDataset, int xsize, int ysize)
{
    GDALRasterBand *poBand1 = poDataset->GetRasterBand(1);
    GDALRasterBand *poBand2 = poDataset->GetRasterBand(2);

    //申请存放NDVI内存
    float *NDVI = (float*) CPLMalloc(sizeof(float)*xsize*ysize);

    unsigned short *pbanddata1 = (unsigned short*) CPLMalloc(sizeof(unsigned short)*xsize);
    unsigned short *pbanddata2 = (unsigned short*) CPLMalloc(sizeof(unsigned short)*xsize);

    //遍历
    for (int i = 0; i < ysize; i++)
    {
        //按行读取图像数据
        poBand1->RasterIO(GF_Read, 0, i, xsize,1,pbanddata1, xsize,1, GDT_UInt16, 0, 0 );
        poBand2->RasterIO(GF_Read, 0, i, xsize,1,pbanddata2, xsize,1, GDT_UInt16, 0, 0 );

        for (int j = 0; j <xsize; j++)
        {
            NDVI[i * xsize + j] = (pbanddata2[j] - pbanddata1[j])/(float)(pbanddata2[j] + pbanddata1[j]);

        }
    }

    //释放内存
    CPLFree(pbanddata1);
    //pbanddata1 = NULL;
    CPLFree(pbanddata2);
    //pbanddata2 = NULL;

    return NDVI;
}

float RSAgricutureDroughtWaterlogging::mean(short *pdata, int nxsize, int nysize)
{
    float avg = 0;
    for (int i = 0; i< nxsize *nysize; i++)
        avg += pdata[i];
    avg /= (nxsize * nysize);
    return avg;
}

float RSAgricutureDroughtWaterlogging::getMax(float *but)
{
    float max_ = but[0];
    for (int i = 1; i < m_ImgWidth * m_ImgHeight; i++)
    {
        if (max_ < but[i])
        {
            max_ = but[i];
        }
    }
    return max_;
}

float RSAgricutureDroughtWaterlogging::getMin(float *buff)
{
    float min_ = buff[0];
    for (int i = 1; i < m_ImgWidth * m_ImgHeight; i++)
    {
        if (min_ > buff[i])
        {
            min_ = buff[i];
        }
    }
    return min_;
}

int RSAgricutureDroughtWaterlogging::saveImageproj(double *adfGeoTransform, const char *proj, unsigned char *buf, const char *outputFileName)
{
    //注册文件格式
    GDALAllRegister();
    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");

    ////注册驱动
    GDALDriver *poDriver;

    ////获取输出文件名后缀，确定输出格式
    std::string outformat = outputFileName;
    char *pszFormat = NULL;
    if (outformat.substr(outformat.size() - 4) == ".tif")
        pszFormat = "GTIFF";//tif格式
    else
    {
        if (outformat.substr(outformat.size() - 4) == ".img")
        {
            pszFormat = "HFA";	//img格式
        }
        else
            return -1;
    }

    //获取支持相应图像的驱动
    poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
    if (NULL == poDriver)
    {
        return -1;
    }

    //创建一个文件，用于保存结果  文件大小为nXsize*nYsize，数据类型为mDataType
    GDALDataset * poDstDS = poDriver->Create( outputFileName, m_ImgWidth, m_ImgHeight,1,GDT_Byte,NULL );

    //设置投影参数
    poDstDS->SetProjection(proj);
    //设置地理转换参数
    poDstDS->SetGeoTransform(adfGeoTransform);


    //保存
    GDALRasterBand  * poBand = poDstDS->GetRasterBand(1);
    poBand->RasterIO( GF_Write, 0, 0, m_ImgWidth, m_ImgHeight, buf, m_ImgWidth, m_ImgHeight, GDT_Byte, 0, 0 );

    //关闭输出文件
    if( NULL != poDstDS )
        GDALClose( (GDALDatasetH) poDstDS );

    //success
    return 0;
}

bool *RSAgricutureDroughtWaterlogging::landsat8Flood(GDALDataset *poDataset, int xsize, int ysize)
{
    //读取波段
    GDALRasterBand *poBand3 = poDataset->GetRasterBand(1);
    GDALRasterBand *poBand4 = poDataset->GetRasterBand(2);
    GDALRasterBand *poBand5 = poDataset->GetRasterBand(3);
    GDALRasterBand *poBand6 = poDataset->GetRasterBand(4);

    //申请存放水域内存
    bool *flood = new bool[xsize * ysize];
    unsigned short *pbanddata3 = (unsigned short*) CPLMalloc(sizeof(unsigned short)*xsize);
    unsigned short *pbanddata4 = (unsigned short*) CPLMalloc(sizeof(unsigned short)*xsize);
    unsigned short *pbanddata5 = (unsigned short*) CPLMalloc(sizeof(unsigned short)*xsize);
    unsigned short *pbanddata6 = (unsigned short*) CPLMalloc(sizeof(unsigned short)*xsize);

    //遍历
    for (int i = 0; i < ysize; i++)
    {
        //按行读取图像数据
        poBand3->RasterIO(GF_Read, 0, i, xsize,1,pbanddata3, xsize,1, GDT_UInt16, 0, 0 );
        poBand4->RasterIO(GF_Read, 0, i, xsize,1,pbanddata4, xsize,1, GDT_UInt16, 0, 0 );
        poBand5->RasterIO(GF_Read, 0, i, xsize,1,pbanddata5, xsize,1, GDT_UInt16, 0, 0 );
        poBand6->RasterIO(GF_Read, 0, i, xsize,1,pbanddata6, xsize,1, GDT_UInt16, 0, 0 );
        for (int j = 0; j <xsize; j++)
        {
            //判断是否为水
            if (pbanddata3[j] + pbanddata4[j] - (pbanddata5[j] + pbanddata6[j]) > 0 &&
                    ((float)pbanddata5[j]) / pbanddata3[j] < 0.88)
                flood[i * xsize + j] = 1;
            else
                flood[i * xsize + j] = 0;

        }
    }
    CPLFree(pbanddata3);
    //pbanddata3 = NULL;
    CPLFree(pbanddata4);
    //pbanddata4 = NULL;
    CPLFree(pbanddata5);
    //pbanddata5 = NULL;
    CPLFree(pbanddata6);
    //pbanddata6 = NULL;

    return flood;
}

bool *RSAgricutureDroughtWaterlogging::modisFlood(GDALDataset *poDataset, int xsize, int ysize)
{
    //读取波段
    GDALRasterBand *poBand1 = poDataset->GetRasterBand(1);
    GDALRasterBand *poBand2 = poDataset->GetRasterBand(2);
    GDALRasterBand *poBand3 = poDataset->GetRasterBand(3);
    GDALRasterBand *poBand4 = poDataset->GetRasterBand(4);
    GDALRasterBand *poBand6 = poDataset->GetRasterBand(5);

    //申请存放水域内存
    bool *flood = new bool[xsize * ysize];
    unsigned short *pbanddata1 = (unsigned short*) CPLMalloc(sizeof(unsigned short)*xsize);
    unsigned short *pbanddata2 = (unsigned short*) CPLMalloc(sizeof(unsigned short)*xsize);
    unsigned short *pbanddata3 = (unsigned short*) CPLMalloc(sizeof(unsigned short)*xsize);
    unsigned short *pbanddata4 = (unsigned short*) CPLMalloc(sizeof(unsigned short)*xsize);
    unsigned short *pbanddata6 = (unsigned short*) CPLMalloc(sizeof(unsigned short)*xsize);

    //遍历
    for (int i = 0; i < ysize; i++)
    {
        //按行读取图像数据
        poBand1->RasterIO(GF_Read, 0, i, xsize,1,pbanddata1, xsize,1, GDT_UInt16, 0, 0 );
        poBand2->RasterIO(GF_Read, 0, i, xsize,1,pbanddata2, xsize,1, GDT_UInt16, 0, 0 );
        poBand3->RasterIO(GF_Read, 0, i, xsize,1,pbanddata3, xsize,1, GDT_UInt16, 0, 0 );
        poBand4->RasterIO(GF_Read, 0, i, xsize,1,pbanddata4, xsize,1, GDT_UInt16, 0, 0 );
        poBand6->RasterIO(GF_Read, 0, i, xsize,1,pbanddata6, xsize,1, GDT_UInt16, 0, 0 );
        //计算均值
        //float bandmean2 = mean(pbanddata2,nXsize,nYsize);
        //float bandmean3 = mean(pbanddata3,nXsize,nYsize);
        for (int j = 0; j <xsize; j++)
        {

            //判断是否为水
            if (pbanddata1[j] - pbanddata2[j] > 0 &&
                    pbanddata1[j] + pbanddata4[j] - (pbanddata2[j] + pbanddata6[j]) > 0 )
                //&&pbanddata2[j] < bandmean2 && pbanddata3[j] < bandmean3)
                flood[i * xsize + j] = 1;
            else
                flood[i * xsize + j] = 0;

        }
    }
    CPLFree(pbanddata1);
    //pbanddata3 = NULL;
    CPLFree(pbanddata2);
    //pbanddata4 = NULL;
    CPLFree(pbanddata3);
    //pbanddata5 = NULL;
    CPLFree(pbanddata4);
    //pbanddata6 = NULL;
    CPLFree(pbanddata6);
    //pbanddata6 = NULL;

    return flood;
}
